Optimal synchronizability of bearings 
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Bearings are mechanical dissipative systems that, when perturbed, relax toward a synchronized 
(bearing) state. Here we find that bearings can be perceived as physical realizations of complex 
networks of oscillators with asymmetrically weighted couplings. Accordingly, these networks can 
exhibit optimal synchronization properties through fine tuning of the local interaction strength as a 
function of node degree [Motter, Zhou, and Kurths, Phys. Rev. E 71, 016116 (2005)]. We show that, 
in analogy, the synchronizability of bearings can be maximized by counterbalancing the number of 
contacts and the inertia of their constituting rotor disks through the mass-radius relation, m ~ r a , 
with an optimal exponent a — a x which converges to unity for a large number of rotors. Under 
this condition, and regardless of the presence of a long-tailed distribution of disk radii composing 
the mechanical system, the average participation per disk is maximized and the energy dissipation 
rate is homogeneously distributed among elementary rotors. 

PACS numbers: 05.45.Xt, 46.55. +d, 45.70.-n, 89.75.-k 



A coherent synchronized motion can naturally emerge 
in a network of oscillators when the coupling intensity ex- 
ceeds the synchronization threshold [1-4 ]. Synchroniza- 
tion is the mechanism responsible for numerous phenom- 
ena such as, e.g., the vital contraction of cells producing 
the heartbeats, the harmony in an orchestra, and the 
coherence of an audience clapping after a performance 
[5H?]. However, undesired synchronization might also be 
responsible for neural diseases and collapse of technical 
infrastructures and networks |S]. Therefore, understand- 
ing how synchronization can be enhanced or mitigated 
is a question of paramount importance. The proper- 
ties of the transition to a synchronized state are known 
to be a result of the interplay between the dynamics of 
the oscillators and the complex topology of the system 
PHU]. Previous studies have shown that synchronization 
can be enhanced on scale-free topologies by asymmetric 
weighted couplings, in contrast to random graphs, where 
the most efficient configuration corresponds to a uniform 
coupling strength [TTJ . More precisely, by expressing the 
interaction strength Si of site i in terms of its degree ki 
as Si = , where (3 is a tunable parameter, Motter 
et at [lTj observed that the properties of the coupling 
Laplacian matrix |12j lead to optimal synchronization at 
P = 1. Under this condition of maximum synchroniz- 
ability, the coupling strength just counterbalances the 
number of connections, thus minimizing the total cost 
associated with the network of couplings. 

Space-filling bearings have been previously considered 
to explain the existence of seismic gaps [2], which are 
those regions between tectonic plates where no earth- 
quake activity has been detected for a large period of 
time [15] . The idea is that the system self-organizes into 
a "bearing state" in which the fragments rotate without 
gliding friction. Systematic procedures have then been 




FIG. 1. Two-dimensional space-filling bearing configuration 
with 31 rotor disks. In order to suppress gliding friction, disks 
rotate either clockwise or anticlockwise, all with the same 
tangential velocity, and loops of touching disks always have 
an even number of disks. This static illusory-motion image 
is an adaptation from the peripheral drift illusion "Rotating 
Snakes" shown in Ref. [13]. Note that every two touching 
disks always have opposite senses of rotation. 



proposed to generate model bearing structures of rotors 
with circular and spherical shapes, in two and three di- 
mensions, respectively, either highly symmetric [161 117] 
or random 18!. As depicted in Fig. [I] hierarchical space- 
filling packings emerge from these models in 2D, where 
the interstices among large disks are sequentially filled by 
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smaller ones. Gliding friction is suppressed by ensuring 
that loops of touching disks have an even number of con- 
stituents. In this way, clockwise turning disks only touch 
counterclockwise rotating ones and vice-versa. At steady 
state, the tangential velocity is the same for all contacts. 

Under a different framework, hierarchically filled struc- 
tures, where smaller elements (e.g, disks and spheres) 
are snugged into the interstices of larger ones, have been 
directly associated to scale-free networks [THHUj- The 
Apollonian packing of circles, for example, inspired the 
introduction of the so-called Apollonian network [T§] . 
where the sites correspond to the centers of the circles, 
and the edges are drawn to connect the centers (sites) of 
pairs of touching circles. Bearings can also be directly 
associated to complex networks, whose sites are given by 
the positions Ri of the centers of the disks. These spatial 
networks fulfill, for each loop of n disks, the condition 

J2"=i (-Sj+i — Ri) = Rn — Ri, and are scale free, if the 
original bearing is space-filling. 

One can readily identify the bearing state as a typical 
synchronized state. It is thus legitimate to convey that 
space-filling bearings rotating in steady-state (i.e., when 
all rotors possess equal tangential velocities) are in fact 
physical realizations of synchronized complex networks. 
Once this conceptual parallelism is ascertained, one can 
go even further and ask, in the spirit of the asymmet- 
ric coupling approach introduced in Ref . [llj , whether or 
not such a synchronized state can be optimized through 
some constitutive physical property of the bearings. In 
what follows we show that the synchronized state of two- 
dimensional space-filling bearings can indeed be substan- 
tially enhanced by adequately adjusting the inertial con- 
tribution of individual rotors to the global motion of the 
system. 

Consider a bearing of N rotors. The equation of mo- 
tion for the angular velocity uji of rotor i can be written 
as, 

h<2% = ^2 Tj = ^2 n^ij x Fji , (1) 
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where h and 7"* are the rotational inertia and radius of 
rotor i, the sum is over all rotors in contact with i, Fji 
is the force of rotor j on the surface of i, and Fjj is the 
unit vector pointing to the contact with j in the reference 
frame of disk i. Taking the force Fji as a dissipative force 
proportional to the relative velocity at the contact point, 
we have, 

Fji = <j(Vj - Vi) = -a(ujj X TjTij + bJ l X TiTij) , (2) 

where a is the coupling between rotors and we used the 
identity fji — —fij. In 2D bearings, since rotors are disks 
with fixed position (see Fig.^J, the net translational force 
is zero and the angular velocity can be described by a 



scalar. Equation (|1J) then simplifies to, 

IiUJi = -a ^ A ij [ u i r f + UJ 3 T i r j\ ) ( 3 ) 
3 

where A^ are the elements of the connectivity matrix, 
defined in such a way that Aij = 1 if two disks i and j are 
different and mutually touching, and Aij — otherwise. 

At this point, we introduce the following constitutive 
relation between the mass of disk i and its radius, 

mi = 2ar? , (4) 

such that the rotational inertia becomes Ii = ar^ +2 , 
where a = 1 in consistent units, for convenience. It then 
follows that, ioi — —P^2jTijU)j, where TJy = riir^ a 5ij + 
r^[ l ~ a rjAij , rii is the number of disks touching i, and 
Sij is the Kronecker delta. In matrix form this can be 
written as, 

'uj = -aTuj , (5) 

where uj is the A-dimcnsional vector of the angular ve- 
locities and T is the interaction matrix. A bearing is 
a dissipative system which, as already mentioned, con- 
verges to a steady-state, namely the bearing state, where 
the tangential velocities of all rotors become equal, i.e., 
Vi = i>2 = • • • = vn = s, such that s(t) = 0. Through 
the relation between tangential and angular velocities, 
uj — R _1 w, where R is a diagonal matrix with Ra = c%ri, 
with Ci = ±1 depending on the sense of rotation of the 
disk [17] , the equivalent to Eq. ^ for the iV-dimensional 
vector of absolute tangential velocities v can be readily 
obtained, being v = — ctBw. The coupling matrix B can 
be written as, B = RTR -1 , i.e., Bij = r~ a (riiSij — Aij), 
where we make use of Ci/cj = — 1 for all pairs of touching 
disks. In this work we focus on the relaxation after small 
perturbations £ to the bearing state, namely Vi = s + £j, 
which leads to the following vectorial variational equa- 
tion: 

i= aBt . (6) 

This system of coupled linear differential equations can 
be written in the space of eigenvectors of B, afjt (with 
eigenvalue Afe), such that X\ (Ai = 0) refers to perturba- 
tions along the stable manifold of the bearing state. All 
other eigenvectors are transverse to x\ [12]. Due to the 
linear nature of Eq. ([6]), the Lyapunov exponents corre- 
spond to — uAfc. Since all eigenvalues are non- negative in 
this problem (the matrix B can be symmetrized) , the sta- 
bility of the bearing state is guaranteed [53] . In analogy 
to the work of Motter et al. [11], the factors r~ a present 
in the elements of B correspond to the weights of the 
pairwise interactions. As we show next, in the complex 
network associated with the topology of the bearings, the 
number of contacts is an increasing function of the radius 
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FIG. 2. Dependence on the exponent a of the lowest (non-zero) A2 (a) and the largest Ajv (b) eingenvalues of the matrix B. 
The same is shown in (c), but for the ratio The insets show the system size dependence for different values of a. Both 

A2 and Ajv increase with a exhibiting a change in the behavior at the crossover value a x (dashed lines): the faster increase 
of A2 is for a < a x , while Ajv grows faster for a > a x . As a consequence, a minimum is observed for the ratio at the 

crossover value a x . 



and so it becomes possible to enhance synchronization by 
tuning the coupling weights in such a way as to balance 
the number of contacts. 

Perturbations along the bearing states (£ = 0) lead 
the system from one bearing state to another with a dif- 
ferent tangential velocity s, i.e., all tangential velocities 
change by the same amount, regardless of the value of a 
(£, = £). These perturbations are related to the eigen- 
value Ai = 0. Hereafter, we focus on perturbations which 
are transverse to the bearing states, i.e., perturbations af- 
ter which the system always relaxes back to the original 
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FIG. 3. a-dependence of the average participation P for 
bearings with different number of disks. An enhanced par- 
ticipation per disk is observed for an intermediate range of 
a. Results suggest a discontinuous change in the position of 
the maximum as discussed in the text. In the inset, the size 
dependence is included for different values of a. The partici- 
pation per disk decreases with the number of disks with the 
smallest slope for a — 1.0. 
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bearing state. Specifically, the objective is to investigate 
the dependence on the inertial parameter a of the small- 
est (non-zero), A2, and largest, A at, eigenvalues of the 
system described by the matrix B. These eigenvalues 
correspond to the slowest and fastest relaxation modes, 
respectively. Generally speaking, for perturbations out 
of the bearing manifold (A ^ 0), the larger the eigen- 
values the faster the relaxation toward the stable state. 
Here two features come into play: the mass (inertia) dis- 
tribution of the disks and the number of contacts. While 
the former depends explicitly on a, the latter is an in- 
creasing function of the rotor radius. Numerical results 
for a bearing (of type n = m = of the first family for 
loops of size 4 [T5]) with 4511821 disks reveal that the 
average number of contacts scales with the disk radius as 
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FIG. 4. Radius dependence of the relative dissipation rate Q 
of each disk for different values of a, obtained for a bearing of 
85159 disks. While for a < 1.0, Q is an increasing function of 
the disk radius r, for a > 1.0 larger disks dissipate less energy 
as the number of contacts does not compensate the increase 
in disk inertia. 
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r 7 , with 7 = 0.94 ± 0.04, which should approach unity in 
the limit of space filling systems. 

As depicted in Figs. [2]ja) and (b), A2 and A^ gener- 
ally increase with a, although changes in the behavior 
of both eigenvalues can be observed at a crossover value 
a x w 1. While An increases faster for a > a x , the in- 
crease of A2 becomes attenuated in the same range of a 
values. The insets of Figs. [2^a) and (b) show that, for 
a > a x , A2 and A^v are increasing functions of the num- 
ber of disks since the ones with higher inertia also possess 
more contacts to shed any perturbation. The same de- 
scription applies for all eigenvalues. For a < a x , the 
inertia assigned to each disk does not always compensate 
its number of contacts. For example, the fact that A2 de- 
creases with N for a = — 1 reflects this imbalance. More 
general conclusions can be drawn from the eigenratio be- 
tween the largest and smallest (non-zero) eigenvalues of 
the coupling matrix, [121 [24]. This ratio solely 

depends on geometrical features (radii, network topol- 
ogy, and disk mass) and not on the initial conditions for 
the velocities. The lower the ratio, the higher the syn- 
chronizability [TH [H] . As shown in Fig. [2^c) , bearings 
consisting of rotors with masses that follow the relation, 
m ~ r Qx , have a minimum eigenratio value, in analogy 
to the optimal synchronization coupling found for scale- 
free networks [TT]. Interestingly, the larger the value of 
N, the more sensitive is the system to variations of a. 
Moreover, the results in the inset indicate that, regard- 
less of the value of a, the eigenratio increases monoton- 
ically with the number of disks. Nevertheless, since this 
increase is less pronounced for a = 1 than for any other 
value of a, the relative depth of the minimum augments 
with system size. As a consequence, large systems display 
an enhanced relative synchronizability at a — a x . In the 
Supplemental Information [25] we show the convergence 
of a x towards unity in the thermodynamic limit. 

To shed light on the dependence on a of the contri- 
bution of disks to the eigenmodes, we define the average 
participation [26] as, 



P = 



TV ^ 



£,-=i**Cj') 4 



(7) 



where the outer sum is over all eigenvectors x k and the 
inner sums are over the components of the eigenvector. 
As shown in Fig. [3j there is an intermediate range of a 
values, for a given system size N, where an enhanced 
average participation per disk can be clearly observed. 
However, our results suggest that the value of a for which 
the participation becomes maximum increases in a dis- 
continuous fashion from a « 0.5 to 1.0, as the system 
size increases from N = 613 to 31531. The inset of Fig. [3] 
shows that, notwithstanding the value of a, the average 
participation per disk decreases with the number of disks 
N, but with a slope that is milder for a — 1.0 than for 
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FIG. 5. Snapshots of bearings with 613 disks and different 
values of a. The color stands for the average dissipation rate 
Q in a logarithmic scale. While for a < 1 the dissipation 
rate is an increasing function of the radius, for a > 1 is a 
decreasing function. For a — 1 the dissipation rate is of the 
same order of magnitude to every disk. 



any other case. This behavior expresses the more homo- 
geneous contribution of rotors to the system dynamics. 
In this range of a values, the number of contacts of each 
disk is approximately compensated by its inertia. Ac- 
cordingly, we expect a x to converge toward 7 rs 1.0 as 
the number of disks increases. Under these conditions, 
the impact of changes in the tangential velocity of a disk 
becomes relatively less dependent on its radius. 

It is also interesting to calculate the rate at which en- 
ergy is dissipated. For rotor i, we have W% = hviVi/r 2 . 
Now, given an eigenmode fc with eigenvalue A&, the ve- 
locity of disk i can be expressed by the i-th component 
of Xk(i) scaled by a constant c, i.e., Vi = cxk{i)- Thus, 
one obtains the relation Vi — —c\k(JXk(i). Based on this 
result, we can decompose the dissipation rate of disk i 
in its different modes, Wki — ~rfc 2 XkCXk(i) 2 . It is then 
possible to quantify how dissipation is distributed among 
disks, by defining the relative dissipation of disk i in the 
eigenmode k as, = Wki/J2j Wkj, which leads to, 



rfx k {i) 2 
E 3 rfx k (j) 2 



(8) 



Figure [4] shows, for different values of a, the dependence 
on the radius r of the relative dissipation rate of a ro- 
tor i averaged over all eigenmodes, Qi = l/AlJ2kQi> 
where M is the total number of eigenmodes. For a < 1.0, 
the dissipation is an increasing function of r, while, for 
a > 1.0, larger disks dissipate less energy than smaller 
ones. By tuning a = 1.0, the dissipation rate becomes 
more uniformly distributed among all disks, being ap- 
proximately invariant on the rotor size. Figure [5] con- 
sists of snapshots for three different values of a, showing 
how the average dissipation rate is typically distributed 
among disks in a bearing with 613 disks. 

In summary, we have shown that bearings are physical 
realizations of complex networks of oscillators. When the 
bearings consist of rotors of different sizes, the coupling 
between oscillators is asymmetric, as the effect of a pair- 
wise interaction on the rotors motion depends on their 
inertia (typically different for each one). Once this par- 
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allelism is established, it is possible to evaluate the stabil- 
ity of these mechanical systems applying concepts from 
dynamic systems theory. In particular, our results for 
two-dimensional space-filling bearings, characterized by 
a scale-free distribution of rotor (disk) contacts, indicate 
that their synchronizability can be duly maximized. This 
is achieved by counterbalancing the number of contacts 
of the disks with their inertia through the mass-radius 
relation, m ~ r"*, where a x is the optimal exponent 
which we expect to asymptotically converge to unity as 
the number of rotors increases. Under this condition, in 
spite of the power-law distribution of radii, the average 
participation per disk has a maximum and the energy dis- 
sipation rate is homogeneously distributed among disks. 
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Fig. SI: Scaling of 1 — a x with the number of disks in 
the bearing, N . A power-law is observed supporting the 
prediction of a x = 1 in the thermodynamic limit. 



